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Abstract 

We consider the problem of stochastic averaging of a quantum two-state dynamics 
driven by non-Markovian, discrete noises of the continuous time random walk type 
(multistate renewal processes). The emphasis is put on the proper averaging over 
the stationary noise realizations corresponding, e.g., to a stationary environment. 
A two state non-Markovian process with an arbitrary non-exponential distribution 
of residence times (RTDs) in its states with a finite mean residence time provides a 
paradigm. For the case of a two-state quantum relaxation caused by such a classical 
stochastic field we obtain the explicit exact, analytical expression for the averaged 
Laplace-transformed relaxation dynamics. In the limit of Markovian noise (implying 
an exponential RTD), all previously known results are recovered. We exemplify new 
more general results for the case of non-Markovian noise with a biexponential RTD. 
The averaged, real-time relaxation dynamics is obtained in this case by numerically 
exact solving of a resulting algebraic polynomial problem. Moreover, the case of 
manifest non-Markovian noise with an infinite range of temporal autocorrelation 
(which in principle is not accessible to any kind of perturbative treatment) is studied, 
both analytically (asymptotic long-time dynamics) and numerically (by a precise 
numerical inversion of the Laplace-transformed averaged quantum relaxation) . 
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1 Introduction 



The influence of a stochastic environment on relaxation and charge transfer 
processes in condensed media [1-3] is a longstanding problem of prime im- 
portance in chemical and statistical physics [4-10]. In this context, exactly 
solvable models are rather rare. The method of a stochastic Hamiltonian by 
Anderson and Kubo [11-13], known also under the label of stochastic Liouville 
equation (SLE) approach [14-22] presents a common methodology that has 
been employed over decades [10,22,23]: It is based on the simplifying repre- 
sentation of a stochastic environment by a classical stochastic field that acts 
on the quantum system of interest without a feedback mechanism. This set 
up thus necessarily inherits some shortcomings. In particular, quantum states 
become asymptotically populated with equal weight as if the environment tem- 
perature were infinite. This shortcoming can be cured in an ad hoc manner by 
adding an extra term to the stochastically averaged Liouville equation to en- 
sure a proper thermal equilibrium. What "proper" means, however, is specific 
to the problem under consideration. Therefore, the mentioned shortcoming 
can be overcome rigorously only within a full quantum-mechanical treatment 
of the total, coupled system-environment dynamics [24,25]. Notwithstanding 
this principal difficulty, the SLE approach remains useful and popular over 
the years. The method is particularly appealing because it allows to extract 
exact results for a number of interesting physical models [17,20,22,23,26-29]. 

In view of the central limit theorem, the use of classical random forces with 
a Gaussian statistics is most frequently employed for this sort of semiclassi- 
cal modeling. The case of Gaussian white noise serves frequently as a simple 
model for the corresponding classical stochastic bath; it formally corresponds 
to a heat bath with an infinite spectrum of excitations. For this noise model 
there exists a number of exact results [17,20,22,23,26,29]. Realistic thermal 
baths possess, however, a spectrum with cut-off at finite energies. This cir- 
cumstance then renders temporal autocorrelations that decay on a finite time 
scale. The use of stationary Gaussian Markov noise with an exponentially de- 
caying temporal autocorrelations (Ornstein-Uhlenbeck process) presents one 
of the simplest models for the corresponding colored noise [20,30]. We remark, 
however, that even for the archetype case of a two-state tunneling system this 
colored noise model cannot be solved exactly, except for some special situa- 
tions, see, e.g., in Ref. [27] for a stochastic Landau-Zener model. In practice, 
one is forced to invoke additional approximations; such as in the case of a 
weakly colored noise the method of cumulant expansions [20,21,31]. It is 
equivalent to an expansion in a small parameter (the Kubo number) which is 
the root mean square (rms) of the fluctuations of the characteristic coupling 
energy (in frequency units) times the autocorrelation time of corresponding 
bath fluctuations [31]. Other perturbation schemes can also be used [30]. In 
particular, the opposite limit of large Kubo numbers (quasi-static fluctuations) 
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also allows for a consistent perturbative treatment [30-32]. 

Nonperturbative approaches are, however, much more appealing. Continuous 
state fluctuations can be approximated with jump-like stochastic processes 
possessing a large number of discrete states. An example presents a discretiza- 
tion of the spatial degree of freedom of a continuous diffusion process. The 
class of Markovian discrete state noises provide a rather general framework 
for a formally exact averaging without using any kind of perturbation expan- 
sion [13,16,33,36-39]. The stochastic path integral approach [34,35], adapted 
for such jump processes [16, 33] is particularly convenient. There exist yet 
other powerful approaches, cf. in the Refs. [36,40-43]. The two-state Marko- 
vian noise (or dichotomous noise) presents the archetype discrete state process 
which allows for an exact study of stochastically driven two-level quantum 
systems [40,44-50]. Furthermore, the multistate case of an exciton transfer 
in molecular aggregates with many independent noise sources modeled by in- 
dependent two-state Markovian noises is also solvable, in the sense that it 
can be reduced to solving numerically an algebraic relation [51]. This dichoto- 
mous noise can serve to model a pseudo-spin 1/2 stochastic bath degree of 
freedom. In the case of electron transfer in molecular systems such pseudo- 
spin stochastic variables can simulate, for example, the bistable fluctuations 
of a charged molecular group nearby the donor, or the acceptor site, or also 
the conformational fluctuations of a bistable molecular bridge [52]. 

The method of a stochastic path averaging can readily be generalized onto non- 
Markovian jump processes of the renewal (or continuous time random walk) 
type [53-55]. Such processes are then characterized by the set of residence 
time distributions (RTDs) ipj{r) in the noise states j and the probabilities pij 
for undergoing a transition from state j — > % at the end of each residence time 
interval (RTI) into another noise state %. The lengths of subsequent RTIs are 
mutually independent. 1 This constitutes the crucial ingredient that allows 
for an exact averaging for such non-Markovian processes [53]. Attention must 
be paid, however, to the stationarity of the process, i.e. to the proper aver- 
aging over the stationary noise realizations in a stationary environment. This 
issue is not trivial and requires special consideration: Such a stationary aver- 
aging is possible only if the average residence times are finite for all discrete 
noise states. The stationary noise averaging can be performed exactly in the 
Laplace-domain for any non-Markovian process of the considered type and 



This does not necessarily imply that the process is Markovian: If the RTD of some 
state is non-exponential the process is clearly non-Markovian and exhibits memory 
effects. This can be understood due to the following reasoning which can be traced 
back to Ref. [56]. For a memoryless, Markovian process the survival probability 
<I>j(r) := ipj(r')dT f in any state j must obviously assume the following property: 
^i( r i + r 2) = &j( T i)^j( T 2)- The only nontrivial solution of this latter functional 
equation reads <F,(t) = exp(— 7j-t), where jj > is the rate to leave the state j. 
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arbitrary quantum dynamics [54,57]. In particular, the Laplace-transform of 
the corresponding quantum propagator can be written in a general analytical 
form [57]. It yields a complex expression which practically can be elaborated 
on explicitly in some special cases only. In this work, we present analytical 
expressions for averaged quantum two-state dynamics driven by a symmetric 
two-state non-Markovian noise with arbitrary RTD possessing a finite mean 
time. Our results generalize the prior results of Refs. [40,44,50] and reduce to 
the latter ones in the particular case of Markovian noise. 

The paper is structured as follows. We first review the general approach and 
present the main result for the Laplace-transformed stationary-averaged quan- 
tum propagator in Sec. 2. A specific application of this general result to the 
Kubo oscillator is given in Sec. 3. Explicit results for a quantum two-state dy- 
namics driven by a two-state non-Markovian noise are provided in Sec. 4. The 
limit of Markovian noise and a special non-Markovian case are studied for the 
two- state quantum dynamics in Sec. 5. The case of manifest non-Markovian 
noise with infinite mean autocorrelation time is considered in Sec. 6. In Sec. 
7 we provide a resume of our results. 



2 Stochastic Liouville-von-Neumann equation and corresponding 
averaged quantum propagator 



To start out, let us consider an arbitrary quantum system with a Hamilton 
operator H[£(t)] which depends on a classical, noisy parameter This 
stochastic process £(£) can take on either continuous or discrete values. Ac- 
cordingly, the Hamiltonian H randomly in time acquires different operator 
values H[£(t)] which generally do not commute, i.e., [#[£(£)], H [£(£')]] ^ 0. 

The prime objective is then to noise-average the corresponding quantum dy- 
namics which is characterized by the Liouville-von-Neumann equation 

= -icmw) , (i) 

for the density operator p{t) over the realizations of the noise £(£). £[£(£)] in 
Eq. (1) denotes the quantum Liouville superoperator, £[£(£)](•) = (•)] 
In other words, the posed challenge consists in evaluating the noise-averaged 
propagator 

rto+t 

(S(t + t,t )) = (TexpH / mr)]dr]), (2) 

J to 

where T denotes the time-ordering operator. 
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2. 1 Non-Markovian vs. Markovian discrete state fluctuations 



We next specify this task for noise assuming N discrete states The noise 
is generally assumed to present a non-Markovian renewal process which is fully 
characterized by the set of transition probability densities ipij{r) for making 
random transitions within the time interval [r, r + air] from the state j to 
the state i. These probability densities are necessarily positive and obey the 
normalization conditions 

£ / Mr)dr=l, (3) 
i=i J0 

for all j = 1,2,..., N. 

The subsequent residence time-intervals between jumps are assumed to be 
mutually uncorrelated. The residence time distribution (RTD) ipj(r) in the 
state j reads 

^(t)=£^(t) = -^. (4) 

i 

The survival probability $j(t) of the state j follows then as 

/oo 
^{r)dr. (5) 

This constitutes the general scheme for continuous time random walk (CTRW) 
theory [58-61]. 

Several particular descriptions used for such non-Markovian processes of the 
renewal type are worth mentioning. The approach in Ref. [53] in terms of 
time-dependent aging rates kij(t) for the transitions from state j to state i 
corresponds to a particular choice, reading 

^•(r) := £^(r)exp[-£ / k i:j (t)dt\. (6) 

i Jo 

The Markovian case corresponds to time-independent transition rates kij(r) = 
const. Any deviation of ^(t) from this Markovian case then in turn yields 
a time- dependence of the transition rates fcjj(r) which amounts to a non- 
Markovian behavior. Furthermore, the survival probability $j(r) in the state 
j is determined by 

N 

^■(r)=exp[-5: / h 3 (t)dt] (7) 

i=l J0 

and Eq. (6) can be recast as 

^■W:=M^;W- (8) 
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The introduction of such time-dependent "aging" rates presents one possibility 
to describe non-Markovian memory effects; it is not unique though. A different 
scheme follows by defining [55]: 



:= Pij(T)ii>j(T) (9) 

with J2iPij( T ) = 1- The physical interpretation is as follows: The process re- 
mains in a state j for a random time interval characterized by the probability 
density At the end of this time interval it jumps into another state % with 

a generally time-dependent conditional probability Pij[j\ Such an interpreta- 
tion is frequently used in the continuous time random walk theory. Evidently, 
any process of the considered type can be interpreted in this way. By equating 
Eq. (8) and Eq. (9) and taking into account that V ; j( r ) := —d&j(r)/dr one 
deduces that the approach in Ref. [53] can be related to that in Ref. [55] with 
the time-dependent transition probabilities 

"^fS) (10) 

and with the non-exponential probability densities tpj(r), 
i.e., Vj(t) = 7 i (r)exp[-/ r 7 i (t)rft] with 7i (r) := E;M r )- 

The description of non-Markovian effects with the time- dependent transition 
probabilities Pij(r), is rather difficult to derive from the sample trajectories 
of an experimentally observed random process £(£). The same holds true for 
the concept of time-dependent rates. These rates cannot be measured directly 
from the set of stochastic sample trajectories. On the contrary, the RTD i^j(r) 
and the time-independent pij (with pa := 0) can routinely be deduced from 
measured sample trajectories, say, in a single-molecular experiment. The study 
of the statistics of the residence time-intervals allows one to obtain the cor- 
responding probability densities ipjij) and, hence, the survival probabilities 
$j(r). Furthermore, the statistics of the transitions from one state into all 
other states allows one to derive the corresponding conditional probabilities 
Pij. From this primary information a complementary interpretation of exper- 
imental data in terms of time-dependent rates fcij(r) can readily be given as 

Mr) = (ID 

if one prefers this particular "language" to describe the non-Markovian ef- 
fects. Moreover, the description with a constant set p^ provides a consistent 
approach to describe stationary realizations of and consequently to find 
the corresponding averaged quantum evolution [57]. It is this reasoning that 
we shall follow in the following. 
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2.2 Averaging the quantum propagator 



The task of performing the noise-averaging of the quantum dynamics in Eq. 
(2) can be solved exactly because we can make use of the piecewise constant 
character of the noise realizations [16,33,53]. Indeed, let us consider the 
time-interval [to,t] and let us take a frozen realization of assuming k 
switching events within this time-interval at the time-instants U, 

to < h < t 2 < ... <t k <t. (12) 

Correspondingly, the noise takes on the values £j , £j k in the time sequel. 

Then, the propagator S(t,t ) reads 

S(t,to) — e - ^^^* - **^ - ^^- 1 ^** - ** -1 ^..^ - ^^"^* 1- * ^ . (13) 

Let us further assume that the process £(£) has been prepared in the state 
jo at t - Then, the corresponding k— times probability density for this noise 
realization reads 

Pk(ij k ,tk]ij k ^tk-i; — ;^i,*i|^o>*o) = ®j k ( t - t k)i J j k j k - 1 (tk-tk-i)--ipj 1 j (t 1 -to) 

(14) 

for k 7^ and Po(£jo> ^o) = *&jo(t ~ to) f° r k = 0. In order to obtain the noise- 
averaged propagator (5 , (t|t , jo)) conditioned on such a nonstationary initial 
noise preparation in the state j one has to average (13) with the probability 
measure in (14) (for k — 0, oo). This task can be readily be performed by use of 
the Laplace-transform method [it will be denoted in the following as A(s) : = 
/ °° exp(— sr)A(r)dr for any time- dependent quantity -A(t)]. The result for 
(S(s\t ,j )) = Io°° exp(-sT)(S(t + r\t ,j ))dr thus reads [53,57] 

(S(s\t ,j )} = j:(Ms)[I-B(s)]~ 1 ) , (15) 

where the matrix-superoperators A(s) and B(s) read in components 

A H (s):=5 kl $,( r ) e -(*+^])^T, (16) 
Jo 

and _ 

_ /*oo 

B k i(s) := / Mr)e~ {s+iCmT dr , (17) 
correspondingly, and I denotes the unit matrix. 

To obtain the stationary noise average of the propagator it necessary to average 
(15) over the stationary initial probabilities pf Q . The averaging over the initial 
distribution alone is, however, not sufficient to arrive at the stationary noise- 
averaging in the case of non-Markovian processes since the noise realizations 
constructed in the way still are generally not stationary. This principal problem 
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is rooted in the following observation: By preparing the quantum system at 
to = in a nonequilibrium state p(0), the noise is taken at random in some 
initial state £ J0 with the probability pf Q (stationary noise). However, every 
time when we repeat the preparation of the quantum system in its initial 
state, the noise will already occupy a (random) state £ J0 for some unknown, 
random time interval t* q (setting a clock at to — sets the initial time for 
the quantum system, but not for the noise, which is assumed to start in the 
infinite past). Therefore, in a stationary setting a proper averaging over this 
unknown time r* o is necessary. The corresponding procedure implies that the 
mean residence time (Tj) is finite, (r,) ^ oo, and yields a different residence 
time distribution for the initial noise state, ijjf\r); namely, it is evaluated 
to read ipj°\r) = $ J (r)/(r ? ) [62]. Only for Markovian processes where $j(r) 
is strictly exponential, does i]/p(r) coincide with ipj(r). Using this ^\t) 
instead of V , i( r ) f° r the first sojourn in the corresponding state and for the 
time-independent set Pij, the noise realizations become stationary [54,57,62]. 
The corresponding expression for the quantum propagator averaged over such 
stationary noise realizations has been derived in Ref. [57], cf. Eqs. (25), (29) 
therein. In a slightly more general form it reads 

(S( 8 )) = (S(8)) M -^(c(8)-A(8)[I-PD(8)]- 1 PA( 8 )) f-, (18) 

where (S(s)) static is the Laplace-transform of the statically averaged Liouville 
propagator 

{S{T)) statlc :=Y,z- imk]T Pt (19) 
k 

pf = lim^oo pj (t) are the stationary probabilities of noise states which are 
determined by a system of linear algebraic equations [54,57], 



V 



St 



(rj) 



Pn 
(Tn) 



(20) 



and P is the matrix of transition probabilities p^, i.e. the "scattering matrix" 
of the random process Furthermore, the auxiliary matrix operators C(s) 
and D(s) in (18) read in components: 

f OO FT 

C kl (s):=5 kl e-(*+^W ^(r')dr'dr (21) 
J o Jo 

and 




We remark here that this very same averaging procedure can be applied to 
any system of linear stochastic differential equations. 
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3 An archetype case: The Kubo oscillator 



A simple but very instructive application of this general procedure is the 
noise-averaging of the Kubo oscillator [13,60]; reading 

x(t) = ie[Z(t)]x(t) . (23) 

This particular equation emerges in various situations such as in the theory 
of optical line shapes, nuclear magnetic resonance [11,13], and also for single 
molecule spectroscopy [63]. In the context of the stochastic theory of spectral 
line shapes [11,13,63], e [£(£)] in Eq. (23) corresponds to a stochastically mod- 
ulated frequency of quantum transitions between the levels of a "two-state 
atom", or transitions between the eigenstates of a spin 1/2 system. 

The spectral line shape is determined via the corresponding stochastically 
averaged propagator of the Kubo oscillator as [13] 

I(uj) = - lim Re\S(-iu + rj)} . (24) 

Note that the limit rj — > +0 in Eq. (24) is necessary for the regularization of 
the corresponding integral in the quasi-static limit (tj) — > oo. Upon identifying 
C[£k] with — tk in Eq. (18) we end up with 



(S(s)} = E - E '"/^T^ (25) 

^ 1 - j)i{s - iej) ( 1 \ 1 - i) n (s - ie n ) p% 



iei \I - PD(s)J im s-ie n (r n ) 

with D nm (s) = S nm ip m (s — ie m ) 2 . The corresponding line shape follows 
from Eq. (25) by virtue of Eq. (24). This result presents a non-Markovian 
generalization of the pioneering result by Kubo [13] for arbitrary iV-state 
discrete Markovian processes. This generalization applies to arbitrary non- 
exponential RTDs i>k{T~), or, equivalently, in accordance with Eq. (11) also for 
time-dependent transition rates kij(r). This generalization was put forward 
originally in Ref. [57] for a particular case, pf = (iA / J2k( T k) , which corre- 
sponds to an ergodic process with uniform mixing (meaning that in a long 
time limit each state j is visited equally often). 



2 Note that the formal solution of another important problem, namely the the 
(first order) relaxation kinetics with a fluctuating rate, pit) = —T[£(t)]p(t) follows 
immediately from (25) upon substitution €j — ► iFj, see in Ref. [64] containing some 
nontrivial non-Markovian examples. 
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Next we apply this result to the case of two-state non-Markovian noise with 
P12 = P21 = 1 and pf 2 = (r li2 )/[(r 1 ) + (r 2 )]. Eq. (25) then yields: 



(£(*))= E \ . , {Tk \ x + (ei_e2) ' 



2 



fc t^ 2 s - ie * ( r i) + ( r 2) (( r i) + ( r 2»( s - *ei) 2 (s - ie 2 ) 

x [i-^i(5-^i)][i-^(s-^ 2 )] ^ (26) 

1-^1(5-^1)^2(5-^2) 

With (26) substituted into (24) one finds the result for the corresponding 
spectral line shape which matches that presented recently in Ref. [65] derived 
therein by use of a different method. Moreover, in the simplest case of Marko- 
vian two-state fluctuations with ^1,2(5) = 1/(1 + ( T i,2)s) and with zero mean, 
(£(0) = ( r i) e i + ( r 2 ) e 2 = 0, this result simplifies further to read 



In (27), a = y(£ 2 (t)) = \€ 2 — e 1 \yJ{T 1 ){T 2 )/({T 1 ) + (r 2 )) denotes the root mean 
squared (rms) amplitude of fluctuations. Moreover, x — ^/2 + ?<Tsinh(6/2) de- 
notes a complex frequency parameter, where v — 1/ (ti) + 1/ (t 2 ) is the inverse 
of the autocorrelation time of the considered process. Its autocorrelation func- 
tion reads (£(£)£(£')) = o" 2 exp(— u\t — t'\). Furthermore, b = ln((r 1 )/(r 2 )) = 
In |e 2 /ei I is an asymmetry parameter. The spectral line shape corresponding 
to (27) has been obtained by Kubo, reading [13,66], 



/M " ^ + ei ) 2 (J+^) 2 + ^- (28) 

Moreover, the expression (27) can readily be inverted into its time domain. 
Note that the corresponding averaged propagator (S(t)) of Kubo oscillator 
[67], i.e., 



(S(t)) 



-xt 



COS 



x 2 t) + 



X 



- x 2 



sm(y/a 2 - x 2 t) 



(29) 



is complex- valued when the process £(£) is asymmetric, i.e. b 7^ 0. This is 
in accordance with the asymmetry of the corresponding spectral line shape, 
I{—oo) 7^ Derived in a different context [68] for the case of a two-state 

Markovian process with a nonvanishing mean and in quite different notation 
an expression equivalent to (29) is known in the theory of single-molecule 
spectroscopy [68-70]. For a symmetric dichotomous process (with 6 = 0) Eq. 
(29) reduces to the expression (6.10) (with ujq = 0) in Ref. [31]. 
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4 Averaged dynamics of a two-level quantum systems driven by 
two-state noise 



Our non-Markovian stochastic theory of quantum relaxation can be further 
exemplified for the relevant case of a two-state quantum system, reading 

H(t) = + E 2 \2}(2\ + ^(0(|1)(2| + |2)(1|), (30) 

which is driven by a two-state non-Markovian stochastic noise £(£) = {A + , A_} 
with corresponding RTDs ip + (r), ip-(r) and the stationary state probabilities 
p± = (t±)/[(t + ) + (t_)]. This noise assumes the normalized stationary auto- 
correlation function 

k(t) := mmmst/(m 2 )st where 6£(t) := £(f) - (£> st is temporal fluctua- 
tion. Its Laplace-transform reads [71-74] 



~k(s) 



1 1 

+ 



i (i -</;_(*)) 



>+> (t^Js 2 l-^ + ( s )^_(s) 
This dichotomic noise possesses the power spectrum Sn(u): 



(31) 



2(A+ - A ) 2 1 
(t+) + (t_) ^ 



1 — ip + {iu))ip-{i(jj) 



(32) 



It causes (dipole) transitions between two states, |1) and |2), and is zero on 
average (a first interpretation). A different interpretation of the considered 
dynamics can also be given when does not vanish on average. Then, for 
A + > A_ > 0, we are dealing with a quantum tunneling dynamics with 
a fluctuating tunneling matrix element, e.g. due to a fluctuating tunneling 
barrier. 

For the considered case of a two state non-Markovian process with pn = P22 = 
0, P12 = P21 = 1 the general result in Eq. (18) can be simplified further. After 
some cumbersome operator algebra we obtain: 



(S(s)) =pfS + (s) + p s ^(s) - {C + (s)+C.(s) 

-[A + (s)B_(s) + A_(s)][I - B + (s)B4s)]-'A + (s) (33) 
-[A_(s)B + (s) + A + (s)][I - B^s)B + (s)]- l A_(s)}, 

where S±(s) denotes the Laplace-transform of the propagator S±(t) = exp(—i£±t) 
with C± := £[A±] corresponding to the fixed noise value £ = A + and £ = A_, 
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correspondingly. Furthermore, C±(s) is given by Eq. (21), C±(s) = C±±(s) 
and 

A±{s) := / $±(r) exp[-(s + *£ ± )r]dr, (34) 

/•oo 

B±(s):= ^j ± (t) exp[-(s + iC ± )r}dr . (35) 

A quantum evolution characterized in Eqs. (33)-(35) by the Liouville operators 
C± does belong to a rather broad class (e.g., for quantum systems with a finite 
number of states) and is not merely restricted to the case of two-state quantum 
dynamics in Eq. (30). 

The archetype model in Eq. (30) does exhibit a rich behavior. In particu- 
lar, it opens a doorway to study the problem of quantum decoherence of a 
two-state atom under the influence of two-state u l/u a " noise that exhibits 
long range time-correlations with a power law decay (for ip(r) possessing a 
long-time algebraic tail, ipij) oc l/r 3 ~ a , < a < 1) [72-74]. Therefore, this 
model constitutes a prominent problem of general interest. Moreover, it relates 
to recent activities that involve decoherence studies for solid state quantum 
computing [75]. It is convenient to express the Hamiltonian (30) in terms of 
Pauli matrices, a z := |1)(1| -|2)(2|, a x := |1)(2| + |2)(1|, a y :=i(|2)<l|-|l)(2|) 
and the unit matrix /, 



H{t) = -He a z + -K(t)a x + -(E 1 + Eh) I, (36) 

where e = (Ei — E 2 )/h. The dynamics of the density matrix of the quantum 
two-state quantum system is then obtained as p(t) — |[7 + J2i= x ,y,z cr i('t)^'i\ 
with components <Ji(t) = Tr(p(t)<3"j). This latter dynamics evolves on a Bloch 
sphere of unit radius (i.e., the corresponding (scaled) magnetic moment is 
conserved, i.e., 3 , \cr(t)\ = 1). Its rate of change obeys, 



&x(t) = -e Q a y (t), 

v y (t)=e a x (t)-£(t)a z (t), (37) 
a z (t)=i(t)a y (t) , 

or a(t) = F[£(i)]a(t) in the vector form, where 



3 This means that each and every stochastic trajectory runs on the Bloch sphere. 
The stochastically averaged Bloch vector dynamics (<?(t)) is, however, "dissipative"; 
i.e. \ {a(t))\ < 1, because (ai(t)) 2 < (af(t)). Thus, the averaged density matrix (p(t)) 
remains always positive in the considered model, independently of the particular 
model used for the stochastic force £(t). 
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-e 
e -£(t) 

£(*) o 



(38) 



The above theory can readily be applied to the noise averaging of a 3-dimensional 
system of linear differential equations (37) over arbitrary stationary realiza- 
tions of with the obvious formal substitution — iC[£(i)] — > F [£(£)]. To- 
wards this goal, we represent the propagators S±(t) = exp(F[A±]t) for the 
fixed static values of noise £ = A + and £ = A_ as matrix expansions over the 
eigen-modes of evolution exp(i\±h), with A±^ = 0, A±^ = A±^ = — 
where Q± = \J^ + A± are the eigen-frequencies of coherent quantum oscilla- 
tions for constant £ = A±. It reads, 



S ± (t) 



E ^exp^t), 

fc=0,l,2 



(39) 



where 



ft 



± 



A 2 ± 



e A± 




[ft 



± 






4 ) 






/ 


4 




eoA± \ 














1 














eoA± 


• A± 


A| 


V 


f2j_ 




W J 



E ^i fc) = '- 



(40) 



fc=0,l,2 



The corresponding Laplace-transformed matrices entering Eq. (33) are repre- 
sented as 



^±( s )= E 

fc=0,l,2 s - zA^ 



a ± ( s )= e i4* 

fe=0,l,2 

>(*) 



fc) l - jj±{s - iX^) 



s - iAi fc) 



s ± ( s )= E i4*ty ± (*-iA£ ; ), 

fc=0,l,2 



(41) 



C ± (a) = E 

fc=0,l,2 



, fe )l -^ ± (s -iAj fc) ) 



(s-iA*?) 2 
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The results in Eqs. (33), (40), (41) form our basis for further studies. 

Let us evaluate it explicitly for the case of symmetric process, ip+ij) = 
ip-(r) = ip(r) : (t+) = (t_) = (r), with zero mean value, A + = — A_ = A. 
We consider the situation where the state "1" is populated initially with the 
probability one: <x 2 (0) = l,a x ^ y (0) = 0. Then, the Laplace-transform of the 
averaged difference of two populations, i.e., (a z (t)) = (pn(t)) — (p22(t)) is de- 
noted as (a z (s)) = (S zz (s)). After some cumbersome manipulations using a 
computer algebra system (MAPLE) to perform multiple matrix operations, 
we end up with the following compact result: 



,~ /xv s 2 + e 2 2A 2 A zz (s) 



where 



A zz (s) = e 2 [l - ^(s)]{(Q 2 - s 2 )(l - 4,{s + Xl)$(s - iSl)) 

- 2iVt s [tp(s + iSl) - $(s - iil)]} 

- A 2 s 2 [l + ^(s)][l - + itt)][l - ^{s - in)], (43) 
B zz (s) =e 2 [l- fa)] [1 + $(s + in)] [1 + $(s - iO)] 

+ A 2 [l + ^(s)](l - rj)(s + iQ)i>(s - iVl)). 

Note that for the considered initial condition, we find (a x (t)) = (c y (t)) = for 
all times because (S xz (s)) = (S yz (s)) = 0. For e = the result in (42)-(43) 
reduces to one for Kubo oscillator (26) with identical ^1,2 (t). Moreover, for 
the Markovian case, tp(s) = 1/(1 + rs), Eq. (42) reduces to 



la is) ) = s 2 + 2us + u 2 + e 2 

KzKn s 3 + 2z/s 2 + (A 2 + e 2 + z/ 2 )s + AV 1 ' 

where v — 2/(r) is the inverse autocorrelation time. This latter result repro- 
duces the result for the averaged populations (pn(s)) = (l/s + (a z (s)))/2 and 
(^22(5)) = (1/s - (a z (s)))/2 in [40,44]. The same result (44) can also be de- 
duced from the known solution for the Markovian case driven by asymmetric 
two-state noise [48] when specified to symmetric noise limit. 

This population difference possesses several remarkable features: The asymp- 
totic difference between populations is zero, (a z (oo)) = ]im s _> (s{a z (s))) = 0. 
In other words, the steady state populations of both energy levels equal 
1/2, independent of the energy difference Hsq. This result can be elucidated 
best in terms of a "temperature" T a of the (pseudo-)spin system. This spin- 
temperature is formally introduced by using for the asymptotic distribution an 
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Ansatz of the Boltzmann-Gibbs form, (p nn (oo)) = exp[—E n /kBT a ]/J2n ex P[~E n /kBT a ] 
Then 4 , 



T CT := ff . .. N (45) 

D \<pn(°o)>/ 

for two-level systems. In accordance with this definition, the result of equal 
asymptotic populations, (£22(00)) = (pn(oo)) = 1/2 can be interpreted in 
terms of an infinite temperature T a = 00. This constitutes a general find- 
ing: a purely stochastic bath corresponds to an apparently infinite temper- 
ature [22,23]. Thus, this stochastic approach modeling the relaxation pro- 
cess in open quantum systems is suitable only for sufficiently high temper- 
atures k B T 3> h\eo\ [22,23]. An asymmetry of unbiased stochastic fluctua- 
tions does not impact this conclusion, see in Ref. [48]. Moreover, the relax- 
ation to the steady state can be either coherent, or incoherent, depending 
on the noise strength and the value of autocorrelation time. In particular, 
an approximate analytical expression for the rate k of incoherent relaxation, 
(pu(t)) = [1 + exp(—kt)]/2, has been obtained in a limit of small Kubo num- 
bers, K := A /v <C 1, which corresponds to a weakly colored noise [31,36]. 
This analytic result reads [40,44,48] 



A 2 u 

* = ^+a- (46) 



The rate exhibits a resonance feature versus v at v = A similar resonance 
feature occurs also in the theory of nuclear magnetic resonance for weakly 
colored Gaussian noise [76]. Note that in Ref. [48] this notable result has been 
obtained for asymmetric fluctuations of the tunneling coupling possessing a 
non- vanishing mean value (£(£)) 7^ 0. This corresponds to a quantum particle 
transfer between two sites of localization which are separated by a fluctuating 
tunneling barrier. A related problem with the inclusion of quantum dissipation 
has been elaborated in [52] within a stochastically driven spin-boson model. 

Yet another interesting solution can be obtained for (a x (s)} with the initial 
condition reading cr x (0) = 1. The Laplace transform of the solution reads 



{x[)) 8(8* + SP) (r)8^ + n^B xx ( S y [ } 

where 

4 This is the standard definition of the temperature of a spin subsystem in nuclear 
magnetic resonance and similar research areas [76]. It is used also to introduce the 
parlance of formally negative temperatures. 
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A xx (s) = [1 - V(s)][l - il>(s + - *p(s - in)], 
B xx (s) = e 2 [l + ^(s)][l - ^(s + - ^{s - tQ)} 

+ A 2 [l - ^(s)](l - ^(s + ifi)^(s - ift)). 



(48) 



The physical relevance of this solution (47) is as follows: In the rotated pseudo- 
spin basis, a x — > a z , b z — > ct x , ct^ — > a^, this problem becomes mathematically 
equivalent to the problem of the derealization of a quantum particle in a sym- 
metric dimer with a tunneling coupling eo under the influence of a dichoto- 
mously fluctuating energy bias £(£). Therefore, it describes the corresponding 
derealization dynamics and, in particular, allows one to determine whether 
this dynamics is coherent or incoherent, depending on the noise characteristics. 

For the Markovian case, Eq. (47) reduces to 5 



Note that the denominators in Eq. (44) and Eq. (49) are different 6 . In the 
more general case of asymmetric Markovian noise, the corresponding denom- 
inator is a polynomial of 6-th order in s, see in [48]. In the considered case 
of symmetric noise it factorizes into the product of two polynomials of 3-rd 
order, namely into those in the denominators of Eq. (44) and Eq. (49). Thus, 
for a general initial condition the relaxation of a two-level quantum system 
exposed to a two-state Markovian field involves six exponential terms. As a 
matter of fact, this seemingly simple, exactly solvable model can exhibit an un- 
expectedly complex behavior even in the Markovian case of a colored two-state 
noise. However, for certain initial conditions, as exemplified above, the general 
solution being a fraction of two polynomials of s simplifies to the results in 
Eq. (44) and Eq. (49). 

In a general case of non-Markovian noise, the analytical solutions in Eqs. (42) 
and (47) can be inverted numerically to the time domain by use of a numer- 

5 The corresponding dynamics also exhibits a resonance feature versus the inverse 
autocorrelation time v within certain limits [50]. 

6 A notable feature is, however, that both corresponding secular cubic equations 
have the same discriminant, D(A,v, eo) = 0, separating the domains of complex 
and real roots. Hence, the transition from a coherent relaxation (i.e. complex roots 
are present) to an incoherent relaxation (i.e. only real roots are obtained) occurs 
at the same values of noise parameters, independently of the initial conditions. 
The corresponding phase diagram separating regimes of coherent and incoherent 
relaxation (judging from the above criterion) has been devised in [50]. It must be 
kept in mind, however, that the weights of the corresponding exponentials are also 
of importance for the characteristic relaxation dynamics. These weights do depend 
on the initial conditions. 




(49) 
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ical Laplace inversion procedure such as the one detailed in Ref. [77] using 
a computer algebra implementation with arbitrary digital precision [78]. A 
quasi-analytical inversion is, however, still possible in specific non-Markovian 
cases. 



5 A simple non-Markovian case 



We consider the following basic case of non-Markovian noise described by a 
bi-exponential RTD 



ipij) = 6a\ exp(— «ir) + (1 — Q)a 2 exp(— air), (50) 

where a lj2 denote two transition rates which can be realized with probabilities 
9 and 1—9, correspondingly. The considered noise possesses the mean residence 
time 



(r) = 0/ai + (1 - 9)/a 2 (51) 
and the mean autocorrelation time r corr := J °° \k(t)\dt reading with k(t) > 0, 

T corr — Cv T corr) (52) 

where C v = ^J(t 2 ) - (r) 2 /(r) is the coefficient of variation of the RTD [73,74] 
and t^} = (t)/2 is the autocorrelation time of Markovian process possessing 
the same mean residence time (t). The ratio T corr /r^} = Cy serves as a 
convenient quantifier for non-Markovian effects. For this basic non-Markovian 
case considered, we find that C v can be large, see in Ref. [64] for details. For 
example, in the limit C, = a i/ a 2 « ^ « 1, C v xs 2/9 3> 1. Such a noise has 
distinct bursting features [64]. 

For this noise, Eq. (42) reduces to a rational function of s, 



- Wy (53) 

where N zz (s) and D zz {s) are specific polynomials of 5-th and 6-th orders, 
correspondingly. The explicit form of these polynomials is not given here since 
it does not provide much physical insight. 

The inversion of (53) to the time domain reads 
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(^W) = E^|exp(r fc t), (54) 

where r k are the complex roots (Re r k < 0) of D zz (s) = and D' zz (s) : = 
dD zz (s)/ds. Even if r k cannot be determined analytically (except for the case 
eo = 0, where the considered rational function can be simplified to the ratio of 
two polynomials of the third and fourth orders, correspondingly), they can be 
found numerically for any set of parameters entering the problem. The prob- 
lem is thereby quasi-analytically solvable. The same is valid for taking more 
terms in the exponential expansion of ip(r). The whole scheme can be easily 
implemented using a computer algebra system like MAPLE or MATHEMAT- 
ICA. In this respect it is pertinent to note that a power law dependence can be 
well approximated by an expansion of the type ip(r) = J2i c i a i ex P( — o^t). In- 
deed, a power law extending over n decades can be satisfactory approximated 
already by an n-term expansion with properly scaled {an} and {q} [79]; see 
also a practical example of a 6-exponential fitting in [80]. Moreover, the nu- 
merical procedures of the Laplace-transform inversion like one in [78] can be 
implemented. An analytical inversion of Eq. (47) to the time-domain has the 
form similar to Eq. (54), but with different polynomials N xx {s) and D xx {s) of 
5-th and 6-th order, correspondingly. 

5.1 Numerical Results 

We consider first the situation with a small Kubo number K = Ar corr <C 1 
where an approximate analytical result, (a z (t)) = exp(-kt), is available with 
the rate k in Eq. (46) [44,48]. This result can also be obtained using the cu- 
mulant expansion method [31]. It does not depend on whether is the stochas- 
tic process under consideration is Markovian or not. This is illustrated with 
Fig. 1. The following parameters are used in calculations presented in Fig. 
1 (units are arbitrary): energy difference between levels eo = 1, noise am- 
plitude A = 0.5. Furthermore, for the Markovian two-state noise we have 
chosen: a± = 100 (here 6 = 1) and for the non-Markovian one: ct\ = 100, 
a 2 = 2000, 9 PS 0.05263157894. Both noises possess the same autocorrelation 
time r corr = 0.005 while the average residence times differ by an order of 
magnitude; the non-Markovian parameter is C v = 10. On the characteristic 
time scale of quantum relaxation the relaxation process in both cases is excel- 
lently described by the approximate analytical result given above. All three 
lines practically coincide in Fig. 1. In this case, the perturbation theory in the 
small Kubo number K works very well. 

It is interesting to note that the initial decay of populations is (beyond the 
discussed exponential approximation) always Gaussian, i.e. \n(a z (t)) oc — t 2 at 
t — > 0. This "Gaussian" regime can be, however, very short and, therefore, it 
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1000 2000 3000 4000 5000 

t 

Fig. 1. Quantum relaxation in the perturbative regime. All three curves practically 
coincide on the considered time scale, see text for details. 

is not readily visible on the characteristic time scale of the relaxation like in 
Fig.l. 

In the opposite limit, K 1, the distinction between the influence of Marko- 
vian and non-Markovian noises possessing the same (mean) autocorrelation 
time becomes rather distinct, cf. Fig. 2. Here, the following noise parameters 
are used: A = 0.5 ct\ = 0.05, 9 = 1 (for the Markovian noise); a,\ = 0.05 
a 2 = 1, ~ 0.05263157894 (for the non-Markovian noise). In both cases, the 
autocorrelation time is the same r corr = 10, while the mean residence times 
are quite different: (r) = 20 for the Markovian case and (r) = 2 for the 
non-Markovian case. 



6 The manifest non-Markovian case 

For a situation in which the autocorrelation time of diverges, the 

corresponding Kubo number K = Ar corr is infinite and any perturbation the- 
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Fig. 2. Quantum relaxation in the non-perturbative regime of a large Kubo number, 
K = 10. The influence of Markovian and non-Markovian noises with the same rms 
and mean autocorrelation time are distinctly different. Coherent features are more 
pronounced for the non-Markovian noise case. Note that in the latter case the overall 
relaxation time is shorter. 

ory is doomed to fail. The theory developed herein, however, allows to obtain 
convergent numerical results by performing the inverse Laplace transforma- 
tion numerically. We illustrate the strength of our theory with the following 
example of manifest non-Markovian noise with an RTD given by 



i>( s ) = --, r~\ — ? — v> ( 55 ) 

l + s{r}g a (sT d ) , 



where 



tanh(W 2 ) 

= —^f2 ( 56 ) 

with < a < 1. Here, (r) in Eq. (55) is the mean residence time and r d is 
a time constant presenting an additional parameter of the distribution. For 
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Td = 0, an exponential distribution is restored. The properties of the RTD 
determined by Eqs. (55), (56) are discussed in detail in [81]. This RTD thus 
possesses three parameters only, but it exhibits an interesting repertoire of 
effects. In particular, it can encompass up to three different interchanging 
power law regimes with the asymptotic power law assuming the form ip{r) oc 

l/ r 2+a T _^ oo [81]. 

The corresponding noise has l/u 1 ~ a feature in its power spectrum at uj — > 
0. Its mean autocorrelation time r corr diverg Therefore, the 

Kubo number K is formally infinite and a corresponding perturbation theory 
is questionable. Our theory delivers, however, exact result for the Laplace- 
transformed noise-averaged quantum relaxation of the excited level popula- 
tion. Moreover, the mean relaxation time can be formally defined as r re i := 
f£°{<7 z (t))dt. With Eq. (42), this time follows as r re/ = \im s ^ {a z (s)) for any 
function ip(s). Moreover, the asymptotic character of relaxation dynamics as 
t — > oo can be found by using the Tauberian theorems of the Laplace-transform 
method and small-s expansion of ip(s) reading (in leading terms) 



i>(s) « 1 - (t)s + As 1+a (57) 

for any ip(r) possessing the finite mean value (r) and the long-time algebraic 
tail ip(r) oc l/r 2+a . For the particular case in Eq. (55), A = (t)t£/3. 

Using (57) in Eq. (42), after some algebra we obtained in the leading order: 



Several remarkable results follow readily. First, (er 2 (oo)) = lim s ^ [s(5" z (s)] = 
0, i.e. the general drawback of the stochastic Liouville approach is preserved. 
Second, the use of a Tauberian theorem [56] in the above equation yields: 



(59) 

For the case in (55), this latter equation is modified as 



~ m^haj (f) ■ (60) 

Remarkably, this result does not depend on the mean residence time (r). The 
tail of the relaxation curve clearly exhibits a power law, (cr z (t)) oc l/t a . 
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The exact result for the Laplace-transformed relaxation in Eq. (42) can be 
reliably inverted numerically due a generalization of the well-known Stehfest 
method [77] in Ref. [78] which requires implementing this method, for example, 
with a computer algebra system like MAPLE (done here) using a sufficiently 
high digital precision. As a "rule of thumb" the number of digits iV used in our 
calculations should correspond to the number of terms taken in the Stehfest 
asymptotical series expansion [77]. iV must be increased until the numerical 
results converge with the required accuracy. For example, to obtain numerical 
data for Fig. 3 we used N = 256. In this figure, the results are numerically 
precise within the corresponding line width. Note that the standard choice 
N — 16 [77] is inadequate to obtain the correct numerical results for the av- 
eraged relaxation in Fig. 3. Following the reasoning in [78] we checked and 
confirmed these numerical considerations for some test functions with known 
results for the "function" - and its known Laplace-transform". Such an im- 
proved Stehfest method presents one of the best numerical Laplace transform 
inversion methods available nowadays (the core of the MAPLE code contains 
just a few lines, by the way). 

For the results in Fig. 3 the following parameters have been used used: e and 
A are the same as in Figs. 1, 2, i.e. e = 1 and A = 0.5; r = 0.01 (like for the 
Markovian case in Fig. 1), a = 0.5 ("1/u; ' 5 " noise) and = 1. It is interesting 
to compare this case with those in Fig. 1 and Fig. 2. Surprisingly, the relax- 
ation dynamics turns out to be initially practically a single exponential with 
small-amplitude quantum coherent oscillations superimposed. About 90% of 
the initial population difference decays exponentially. The long time tail of 
the relaxation process is, however, clearly non-exponential, and the approach 
to the steady state occurs much slowly than the initial exponential decay. The 
emergence of such a slow non-exponential asymptotic decay is rather intrigu- 
ing. It is due to a manifest non-Markovian character of the noise as detailed 
analytically above. The same reasoning holds true is valid for a = 0.1; the 
corresponding noise obeys an l/u 0SJ feature in its power spectrum which is 
close to 1/f noise. In this case, the tail of relaxation curve becomes, however, 
more flat, being in accordance with Eq. (60), and the approach to the steady 
state (cr z (oo)) = occurs extremely slow: This might create an incorrect im- 
pression that (a z (oo)) ^ 0, see in Fig. 4. We remark that for this particular 
case the digital precision N = 32 was sufficient to obtain convergent results 
within the improved Stehfest method. 



7 Resume 

With this work we have presented general results for the averaged quantum 
relaxation that is driven by discrete state non-Markovian noise of the renewal, 
or CTRW type. The noise sources present non-Markovian generalization of dis- 
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Fig. 3. (a) Quantum relaxation under the influence of non-Markovian two state noise 
with extreme long time correlations. The parameters of this manifest non-Markovian 
two-state noise are: A = 0.5, a = 0.5, (r) = 0.01, Td = 1. (b) The same data 
presented on log-linear scales. 
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Fig. 4. Quantum relaxation under the influence of non-Markovian two-state noise 
with a l/o; 0,9 feature for its power spectrum, a = 0.1, the remaining parameters are 
the same as in Fig. 3. 

crete state Markovian noise sources. Our focus has been on the averaged time- 
evolution in presence of stationary noise realizations (i.e. the noise does not 
relax to, but is in its steady state while the noise-averaged quantum dynam- 
ics undergoes a relaxation process). The practical feasibility of our approach 
has been elucidated with several applications of general interest, namely (i) 
the averaging of the Kubo oscillator (for arbitrary noises) and (ii) the aver- 
aged relaxation dynamics of a quantum two state system that is driven by 
two state non-Markovian noise. For this latter case and for a symmetric pro- 
cess, tractable analytical expressions have been obtained in Eqs. (42), (47). 
The previously known results for the case of Markovian noise are recovered 
from our more general expressions. These new analytical results have been 
verified by corresponding numerical studies. In particular, the case of a mani- 
fest non-Markovian noise possessing an infinite variance of the residence time 
distribution, an infinite mean autocorrelation time and an l/u a feature in 
its power spectrum has been studied. This latter situation cannot be tackled 
within perturbation theory. However, our theory allows for a definite non- 
perturbative treatment. The authors share the confident belief that this new 
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theory will prove useful for this and similar future investigations of quantum 
relaxation processes that are exposed to one or several noisy environments 
exhibiting characteristic long-time memory. 
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